As is known, determination of the best spacecraft trajectory for a given mission is defined in literature as the “Two Point Boundary Problem”, and it consists on a determination of a trajectory, among all possible ones, connecting two different points in space (representing the boundary conditions) and allowing to maximize or minimize a given cost function (also named cost index). Equations of motion are differential constraints of the problem.
Several optimization methods are known to be useful to solve optimization problems related to space transfers, which fall essentially into two main categories:
classical, indirect methods based on Pontryagin maximum principle; and
direct methods, which attempt to find the minimum of the cost function through many numerical procedures and can be further grouped into                direct methods applied to the calculus of variation, and        search techniques.        
Also genetic algorithms are employed to solve problems related to interplanetary transfers.
A brief description of the previously listed known methods is given below.
The Pontryagin maximum principle is a basic theorem in the calculus of variation applied to optimal control theory: it gives some necessary conditions to determine an optimal solution to the problem under analysis, and it is based solely on differential properties that certain classes of functions exhibit at points of extrema. According to such conditions, it is possible to determine an evolution in time of some parameters, the so-called Lagrange multipliers, used to evaluate the control variables needed to solve the optimal problem. The initial values of the Lagrange multipliers are unknown, and they must be evaluated numerically so that boundary conditions are satisfied. In other words, the optimization problem is completely transferred to the determination of the starting values of the Lagrange multipliers. The problem becomes even more complex when the final state is not fully determined, that is, some variables of the state vector are not assigned at the boundary conditions. Beside the conditions given by the Pontryagin maximum principle, additional constraints defined as “transversality conditions” must be considered.
In detail, the indirect methods, through the conditions stated by the Pontryagin maximum principle, allow to determine at runtime which is the best way to employ the control variables to solve the optimization problem, but such variables are unknown at the beginning, as shown below.
Let{dot over (x)}=ƒ(x,u,t)be a system of partial differential equation of the first order, in which x is the state vector and u the so-called control vector, that is the vector to be determined at runtime to minimize or maximize a given cost function J. By recalling the formulation of the Bolza problem, the cost function may have the form:
  J  =                    ∫                  t          0                          t          f                    ⁢                        L          ⁡                      (                          x              ,              u              ,              t                        )                          ⁢                                  ⁢                  ⅆ          t                      +          g      ⁡              [                              t            0                    ,                      x            ⁡                          (                              t                0                            )                                ,                      u            ⁡                          (                              t                0                            )                                ,                      t            f                    ,                      x            ⁡                          (                              t                f                            )                                ,                      u            ⁡                          (                              t                f                            )                                      ]            
Let us assume, for example, the cost function J must be minimized. This leads to the following condition:
      J    *    =            inf              {                  u          ∈          U                }              ⁢    J  
Let also H(x,u,λ,t) be the Hamiltonian function associated to the system, written in the form:H(x,u,λ,t)=λT{dot over (x)}−L. 
Note that if another given functional l, defined as
  I  =                    ∫                  t          0                          t          f                    ⁢                        M          ⁡                      (                          x              ,              u              ,              t                        )                          ⁢                                  ⁢                  ⅆ          t                      +          h      ⁡              [                              t            0                    ,                      x            ⁡                          (                              t                0                            )                                ,                      u            ⁡                          (                              t                0                            )                                ,                      t            f                    ,                      x            ⁡                          (                              t                f                            )                                ,                      u            ⁡                          (                              t                f                            )                                      ]            has to be maximized, the previous formulas would become:
      I    *    =            sup              {                  u          ∈          U                }              ⁢    I  andH(x,u,λ,t)=λT{dot over (x)}+M, where λ is the vector of Lagrange multipliers (which is also defined as the co-state vector or the adjoint variable vector). The Pontryagin maximum principle states that if x(t), λ(t) satisfy the conditions:
            x      .        =                  ∂        H                    ∂        λ                        λ      .        =          -                        ∂          H                          ∂          x                    and for all tε└t0;tƒ┘H(x,u,λ,t)≦H(x,u•,λ,t),then u• is the desired optimal control which maximize the Hamiltonian function H.
Beside the Pontryagin maximum principle, if the state vector is not fully determined at the boundaries, farther conditions named “transversality conditions” must be added. The number of these additional constraints equals the amount of entries in the state vector which are undetermined at the boundaries. For example, if the final time is not established, the corresponding transversality condition shall be:Hƒ=1if the minimum transfer time problem is considered, or:Hƒ=0if the minimum propellant consumption problem is under analysis.
The number and the expression of the transversality conditions changes according to the given problem under analysis.
Let us apply the general formulation above to the two points problem, in the case of interplanetary transfers, for example let us consider an Earth-Mars transfer. Under the hypothesis of a co-planar transfer (all the orbits of the major bodies in the Solar System have a small inclination with respect to the Ecliptical plane), and considering an inertial reference frame xOy with the origin O coincident with the center of mass of the Sun, the motion equations can be written as follows:
         {                                                    ⁢                                                            ⅆ                  x                                                  ⅆ                  t                                            =              u                                                                                    ⁢                                                            ⅆ                  y                                                  ⅆ                  t                                            =              v                                                                                    ⁢                                                            ⅆ                  u                                                  ⅆ                  t                                            =                                                                    -                                          μ                                                                        (                                                                                    x                              2                                                        +                                                          y                              2                                                                                )                                                                          3                          2                                                                                                      ⁢                  x                                +                                  δ                  ⁢                                      T                    m                                    ⁢                                      v                    x                                                                                                                                    ⁢                                                            ⅆ                  v                                                  ⅆ                  t                                            =                                                                    -                                          μ                                                                        (                                                                                    x                              2                                                        +                                                          y                              2                                                                                )                                                                          3                          2                                                                                                      ⁢                  y                                +                                  δ                  ⁢                                      T                    m                                    ⁢                                      v                    y                                                                                                                                    ⁢                                                            ⅆ                  m                                                  ⅆ                  t                                            =                                                -                  δ                                ⁢                                                    k                                                                                          where {x,y,u,v,m} is the state vector (x,y are the position vector components whereas u,v are those of the velocity vector; m is the actual mass of the spacecraft), μ is the Sun gravitation parameter, T a thrust level achievable with a given propulsion system (continuous manoeuvre) and δ is an on-off function representing the propulsion system state (δ=0, engine is off; δ=1 engine is active). {vx,vy} is the unit vector indicating the thrust direction. Given a departure time to and an arrival time tf of the spacecraft, let also be the functional
  J  =                    t        f            -              t        0              =                  ∫                  t          0                          t          f                    ⁢              1        ⁢                                  ⁢                  ⅆ                                          ⁢          t                    that has to be minimized (i.e., transfer time has to be minimized); the corresponding Hamiltonian function shall assume the form:
  H  =                    λ        x            ⁢      u        +                  λ        y            ⁢      v        +                  λ        u            ⁡              [                                            -                              μ                                                      (                                                                  x                        2                                            +                                              y                        2                                                              )                                                        3                    2                                                                        ⁢            x                    +                      δ            ⁢                          T              m                        ⁢                          v              x                                      ]              +                  λ        v            ⁡              [                                            -                              μ                                                      (                                                                  x                        2                                            +                                              y                        2                                                              )                                                        3                    2                                                                        ⁢            y                    +                      δ            ⁢                          T              m                        ⁢                          v              y                                      ]              -                  λ        m            ⁢      δ      ⁢                      k                      -    1  or in a more useful way:
  H  =                    λ        x            ⁢      u        +                  λ        y            ⁢      v        +                  λ        u            ⁡              [                              -                          μ                                                (                                                            x                      2                                        +                                          y                      2                                                        )                                                  3                  2                                                              ⁢          x                ]              +                  λ        v            ⁡              [                              -                          μ                                                (                                                            x                      2                                        +                                          y                      2                                                        )                                                  3                  2                                                              ⁢          y                ]              +          δ      ⁡              [                                            T              m                        ⁢                          (                                                                    λ                    u                                    ⁢                                      v                    x                                                  +                                                      λ                    v                                    ⁢                                      v                    y                                                              )                                -                                  k                                      ]              -    1  in which {λx, λy, λu, λv, λm} is the Lagrange multiplier vector (co-state vector). Let's apply the Pontryagin maximum principle. The control is optimal if and only if:
         {                                                      v              x                        =                                          λ                u                                                                                  λ                    u                    2                                    +                                      λ                    v                    2                                                                                                                                      v              y                        =                                          λ                v                                                                                  λ                    u                    2                                    +                                      λ                    v                    2                                                                                          and
  δ  =            1      ⁢                          ⁢              if        ⁢                                  [                                            T              m                        ⁢                          (                                                                    λ                    u                                    ⁢                                      v                    x                                                  +                                                      λ                    v                                    ⁢                                      v                    y                                                              )                                -                                  k                                      ]              >    0  
The first condition states that the thrust direction vector {vx,vy} must be aligned with the {λu,λv} vector, whereas the second condition gives an indication whether it is convenient to set the thruster on or not. As for the co-state derivatives, the Pontryagin maximum principle states that:
         {                                                    ⁢                                                            ⅆ                                      λ                    x                                                                    ⅆ                  t                                            =                                                μ                                                            (                                                                        x                          2                                                +                                                  y                          2                                                                    )                                                              5                      2                                                                      ⁡                                  [                                                                                    λ                        u                                            ⁡                                              (                                                                              2                            ⁢                                                          x                              2                                                                                -                                                      y                            2                                                                          )                                                              +                                                                  λ                        v                                            ⁡                                              (                                                  3                          ⁢                          xy                                                )                                                                              ]                                                                                                                  ⁢                                                            ⅆ                                      λ                    y                                                                    ⅆ                  t                                            =                                                μ                                                            (                                                                        x                          2                                                +                                                  y                          2                                                                    )                                                              5                      2                                                                      ⁡                                  [                                                                                    λ                        v                                            ⁡                                              (                                                                              2                            ⁢                                                          y                              2                                                                                -                                                      x                            2                                                                          )                                                              +                                                                  λ                        u                                            ⁡                                              (                                                  3                          ⁢                          xy                                                )                                                                              ]                                                                                                                  ⁢                                                            ⅆ                                      λ                    u                                                                    ⅆ                  t                                            =                              -                                  λ                  x                                                                                                                  ⁢                                                            ⅆ                                      λ                    v                                                                    ⅆ                  t                                            =                              -                                  λ                  y                                                                                                                  ⁢                                                            ⅆ                                      λ                    m                                                                    ⅆ                  t                                            =                                                T                                      m                    2                                                  ⁢                                  (                                                                                    λ                        u                                            ⁢                                              v                        x                                                              +                                                                  λ                        v                                            ⁢                                              v                        y                                                                              )                                                                        
If all these conditions are respected, the Pontryagin maximum principle is fully satisfied, and the optimal thrust direction is given analytically through the co-state variables λu,λv. The main problem affecting this approach is that the initial values of {λx, λy, λu, λv, λm} is unknown and it cannot be estimated through physical consideration; it is only possible to suppose the order of magnitude of each term by looking at the corresponding derivative. By far, there is a limited number of optimal vectors {λx,0, λy,0, λu,0, λv,0, λm,0} that must be searched in  (or in  if the hypothesis of coplanar transfers lapses) allowing to reach the desired final conditions. It is clear that such an approach, although it gives an analytical expression for the optimal thrust direction and magnitude, appears quite unsatisfying because of the difficulty to reach the final conditions.
Instead, the direct methods applied to the calculus of variation are based on the approximation of the state vector and of the cost-functional through series representations. Once a proper set of elementary functions has been chosen for the series representation, under given hypotheses the unknown coefficients of the series can be found by imposing the extremization of the approximate cost functional. In other terms, the problem of extremizing a functional with respect to a set of functions is reduced to the extremization of a function with respect to a set of parameters. The solution obtained is in general an approximate solution, unless the set of parameters is an infinite set.
In detail, this kind of direct methods aims to reduce the problem of maximizing (or minimizing) a functional with respect to a function to the problem of maximizing (or minimizing) a function with respect to a set of variables, through series representations, albeit, the solution obtained is in general an approximate solution, unless the set of parameters is an infinite set. Let's consider a 1-dimensional differential problem written in the form {dot over (x)}=ƒ(x,u,t), in which x is the state vector and u the control vector. Let also J(x) be the functional to be maximized or minimized. Series approximation methods are based on the assumption that a function x=x(t) which makes maximum (or minimum) a functional J(x) can be expanded in a series of the form:
      x    ⁡          (      t      )        =            ∑              k        =        0            ∞        ⁢                  a        k            ⁢                        ψ          k                ⁡                  (          t          )                    and also
            x      .        ⁡          (      t      )        =            ∑              k        =        0            ∞        ⁢                  a        k            ⁢                                    ψ            .                    k                ⁡                  (          t          )                    where ak's are coefficients independent of time and Ψk(t) are known functions. It is clear that truncated forms of the previous expressions must be employed:
                    x        n            ⁡              (        t        )              =                  ∑                  k          =          0                n            ⁢                        a          k                ⁢                              ψ            k                    ⁡                      (            t            )                                                            x          .                n            ⁡              (        t        )              =                  ∑                  k          =          0                n            ⁢                        a          k                ⁢                                                            ψ                .                            k                        ⁡                          (              t              )                                .                    
It is assumed that this sequence of approximating functions is complete for functions x(t). By definition, a sequence of approximating functions is complete under the following condition: for all ε>0, there exists an integer j such that:
            min                        a                      k            ′                          ⁢        s              ⁢          [                        max                      t            ∈                          [                                                t                  0                                ;                                  t                  f                                            ]                                      ⁢                                                                      x                j                            ⁡                              (                t                )                                      -                          x              ⁡                              (                t                )                                                                  ]        <  ɛThe importance of this definition lies in the fact that
            min      x        ⁢                  ⁢          J      ⁡              (        x        )              =            min                        a                      k            ′                          ⁢        s              ⁢          [                        lim                      j            →            ∞                          ⁢                                  ⁢                  J          ⁡                      (                          x              j                        )                              ]      if J(x) is strongly continuous and if the sequence of functions is complete. The approximate functional Jn(xn) shall assume the form:
            J      n        ⁡          (              x        n            )        =            ∫              t        0                    t        f              ⁢                  L        ⁡                  (                                    x              n                        ,                                          x                .                            n                        ,            t                    )                    ⁢              ⅆ        t            
To solve the optimization problem, Jn(xn) must be stationary with respect to each ak's, that is:
            ∂                        J          n                ⁡                  (                      x            n                    )                            ∂              a        k              =            ∫              t        0                    t        f              ⁢                            ∂                      L            ⁡                          (                                                x                  n                                ,                                                      x                    .                                    n                                ,                t                            )                                                ∂                      a            k                              ⁢              ⅆ        t            
Once, the indicated integrations are effected, the resulting set of equations can be solved for stationary points, and those stationary points which are relative maxima (or minima) can be determined through search techniques. The functional Jn(xn) is often expressed as follows:
  f  =            ∑              j        =        0            m        ⁢                  ∑                  i          =          0                m            ⁢                        b          ij                ⁢                  x          in                ⁢                              x            .                    jn                    
Moreover, the search techniques are based on numerical procedures to evaluate the optimal solution starting from the approximate evaluation of the cost function. Once such an approximation is obtained, the set of variables maximizing (or minimizing), the cost function are then searched in a given domain; therefore such methods are also defined “search techniques”. The function can be given analytically or it can be determined experimentally; it may or may not exhibit discontinuities, and there can be constraint equations limiting the arguments of the performance measure. In the latter case, the problem is called “the non linear programming” problem. Among the maximum (or minimum) search techniques, gradient-based methods, univariate search and non-sequential search are worth to be mentioned.
Furthermore, the genetic algorithms are search techniques used to find exact or approximate solutions to optimization and search problems. They belong to the wide class of evolutionary algorithms, that use theoretical constructs inspired by evolutionary biology such as inheritance, mutation, selection, and crossover. Several sets consisting of two or more control parameters, called chromosomes, are chosen to solve the optimization problem; once the chromosomes have been built, a random “matches” between chromosomes take place, and through the crossover, new chromosomes (offspring) hopefully closer to the desired solution are breed. Crossovers occurs according to given rules. The whole process is then repeated until no significant improvements are obtained; in this way the algorithm finds the optimal solution.
All the previously described optimization techniques are affected by intrinsic limitations, bounding their pertinence to a given problem or their efficiency in finding a solution.
Furthermore, the Applicant has noted that the “Two Point Boundary Problem” applied to spacecraft transfers, not only interplanetary ones, has proved to be more complex than it was thought to be (this result will be described below), and this strongly affect the validity of each previously described method.
As for the indirect methods, although some analytic constraints for optimality are satisfied, there is no way to establish a priori the value of the initial Lagrange multipliers allowing to meet the boundary conditions. A three dimensional analysis of the optimal two point problem would require six multipliers for the spacecraft state vector, plus one corresponding to its mass. Furthermore, if the initial epoch is a parameter for optimization, also the starting instant t0 to begin the mission must be considered. About the initial epoch t0, it may be important to underline that, when considering interplanetary transfers, especially when transfers among inner planets are under analysis, the choice of the initial epoch to plays a fundamental role in the mission planning, since the epoch rules the mutual position between the starting point and the target. If the mutual position is not chosen properly, the space mission might also fail.
Therefore, as for the indirect methods, up to eight variables for optimization can be included in the analysis for optimization, and this leads to a large expense of computational time.
As for the direct methods applied to the calculus of variations, it is quite difficult to obtain an analytic expression for the error induced by the previously described procedure, therefore it is difficult to estimate whether the approximation of the cost function converges or not, so it is difficult to understand if the approximate solution matches the real behavior of the optimal solution. Furthermore, the analysis of N-dimensional problems, such as an interplanetary transfer, requires a wide number of integrations and coefficients, but all these efforts might not be sufficient by themselves to ensure a satisfying representation of the cost functional because of the cross-dependence of each variable in the state vector.
Another limitation affecting this kind of methods is that the cost functional, as well as the state vector components, must be continuous in the domain of integration. Problems involving instantaneous changes in the state vector, for example the switch of reference frame during an interplanetary transfer or the loss of a stage when analyzing a launcher trajectory, cannot be approached reliably through these methods.
The gradient methods are based on the determination of the cost function gradient through several evaluations of the function itself with respect to given parameters; the better is the knowledge of the function gradient, the greater is the probability of finding a minimum (or maximum) in the search domain. Local minimums can be obtained quickly if the cost function is smooth with respect to the chosen parameters; yet problems such as spacecraft trajectory evaluation are characterized by cost functional which is striking irregular, as shown in the followings. Therefore gradient-based search techniques might not be reliable when considering this kind of problems.
The univariate search methods are quite simply and easy to be implemented. Such methods have been used, for example, by everyone who has had occasion to tune an electrical circuit by adjusting several parameters. First, one parameters is adjusted until no farther improvement is gained; then another parameter is tuned until no additional improvement results, and so on. After each parameter has been adjusted once, the process is repeated by returning to the first parameter and proceedings as described before. Under the non-plausible hypothesis that the parameters do not interacts one with each other, this procedure leads to the desired optimal configuration. Cross-dependence among variables strongly affects fields of application of this procedure.
The non-sequential search is conducted over evenly spaced points in a simply connected region of a Euclidean space. Each of the xi coordinates is assigned a set of evenly spaced points, called grid points, and only the values of xi at these grid points are used. The cost function is then evaluated for all the possible combinations of the grid points, and the grid value of the vector x leading to the best ƒ(x) is deemed the winner. The grid spacing must be chosen accordingly by the searcher; as for the number of data points required, if each xi coordinate is assigned k spaced points, the total number of data is kn, being n the number of entries in x. The non-sequential search techniques do not require assumptions concerning the abruptness of changes in the state vector and in the cost functional, but the large number of evaluations required makes them quite inefficient.
Let's, now, consider the Two Point Boundary Problem, and assume that the overall transfer time is the functional to be minimized. A general technique used to solve the two points problem is to minimize an
            J      ~        ⁡          (                        t          f                ,                  r          target                ,                  V          target                ,                  γ          f                    )        =                                          [                                          (                                                      r                    f                                    -                                      r                    target                                                  )                                            δ                ⁢                                                                  ⁢                r                                      ]                    2                +                              [                                          (                                                      V                    f                                    -                                      V                    target                                                  )                                            δ                ⁢                                                                  ⁢                V                                      ]                    2                +                              [                                          (                                  γ                  f                                )                                            δ                ⁢                                                                  ⁢                γ                                      ]                    2                      +          [                        t          f                          t          ref                    ]      augmented cost function {tilde over (J)}, which is a weighted combination of the cost functional itself plus the norm of the error evaluated at the boundaries:where (rf−rtarget) is a difference between the spacecraft final position vector and the target position vector, (Vf−Vtarget) is a difference between the spacecraft final velocity vector and the target velocity vector, and γƒ is the flight angle, which must be close to zero so that the spacecraft transfer trajectory is eventually tangent to the target orbit at the encounter. Coefficients δr,δV,δγ represent tolerances, as well as tref, that have been introduced to make the functional homogeneous. Such a cost functional is useful to increase the likelihood of convergence towards the optimal solution with the given boundary constraints. The problem arising is that this kind of function shows two distinct trends which make finding the optimal solution very difficult.
In particular, this kind of function shows a smooth, macroscopic trend and an irregular, microscopic behaviour, as shown in FIG. 1 and in FIG. 2, respectively.
In detail, FIG. 1 shows an example of the macroscopic trend of the augmented cost function {tilde over (J)} which appears smooth in wide regions of space, while FIG. 2 shows an example of the microscopic trend of the augmented cost function {tilde over (J)} wherein the region close to the origin is crowded with a striking number of peaks and valleys.
It represents a proof of the fact that the Two Point Boundary Problem is often much more difficult than it was thought, and the convergence to the true optimal solution is an extremely improbable event, at least through the known optimization algorithms described above.